#Importing packages
#Numpy for vector calculations and plotting
import numpy as np
#Matplotlib for plotting (and animating)
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
#Sympy for vector calculations
import sympy as sym
sym.init_printing()
#Lambdify for converting sympy functions/values to numpy ones
from sympy.utilities.lambdify import lambdify
Written Work:
After determining the sum of forces acting on each mass (as seen above), we get the following matrices for $\textbf{K}$ and $\textbf{M}$:
#Defining matrices K and M
#Defining symbols
k, m = sym.symbols('k, m')
#Defining matrix M
M = sym.Matrix([[m, 0, 0, 0],[0, m, 0, 0],[0, 0, m, 0],[0, 0, 0, m]])
M #displaying result
The above result is our $\textbf{M}$ matrix for our system
#Defining matrix K
K = sym.Matrix([[2*k, -k, 0, 0],[-k, 2*k, -k, 0],[0, -k, 2*k, -k],[0, 0, -k, 2*k]])
K #displaying result
The above result is our $\textbf{K}$ matrix for our system
We can then calculate the normal frequencies through taking the determinant of $\textbf{K}-\omega^2\textbf{M}$, setting the result equal to zero, and solving for $\omega^2$:
#Calculating normal frequencies
#Defining new symbol
ω = sym.symbols('ω')
#Defining our relationship
KwM = K - ω**2*M
KwM #Checking result
#Calculating determinant of above matrix
det = sym.det(KwM)
total_norm_freq = sym.solve(det, ω**2) #Solving for ω^2
#Saving these determinants as variables
ω_norm1 = total_norm_freq[0]
ω_norm2 = total_norm_freq[1]
ω_norm3 = total_norm_freq[2]
ω_norm4 = total_norm_freq[3]
total_norm_freq #Displaying result
Thus our normal frequencies are: $$\omega^2_1 = \frac{k(3/2-\sqrt{5}/2)}{m}$$ $$\omega^2_2 = \frac{k(5/2-\sqrt{5}/2)}{m}$$ $$\omega^2_3 = \frac{k(5/2+\sqrt{5}/2)}{m}$$ $$\omega^2_4 = \frac{k(3/2+\sqrt{5}/2)}{m}$$
Plugging in each of these normal frequencies to the original equation $\textbf{K}-\omega^2\textbf{M}$, as outlined in Part b, we can calculate their corresponding eigenvectors to gather their corresponding normal mode:
#Defining new symbols
a1, a2, a3, a4 = sym.symbols('a_1, a_2, a_3, a_4')
#Defining general matrix - used in determining
A_matrix = sym.Matrix([a1,a2,a3,a4])
#First Normal Mode - substituting first solution from Part B into matrix, multiplying by general matrix,
#and getting relationships between variables a1, a2, a3, and a4 based on resulting system of equations
arg1 = sym.solve(KwM.subs(ω**2,ω_norm1)*A_matrix,(a1,a2,a3,a4))
mode_norm1 = A_matrix.subs(arg1).subs(a4,1) #substituting above arguments into general matrix and evaluating when a4 = 1
mode_norm1 #Displaying result of first normal mode
#Second Normal Mode - same process as above
arg2 = sym.solve(KwM.subs(ω**2,ω_norm2)*A_matrix,(a1,a2,a3,a4)) #same method of generating arguments
mode_norm2 = A_matrix.subs(arg2).subs(a4,1) #same substitution
#Third Normal Mode - same process as above
arg3 = sym.solve(KwM.subs(ω**2,ω_norm3)*A_matrix,(a1,a2,a3,a4)) #same method of generating arguments
mode_norm3 = A_matrix.subs(arg3).subs(a4,1) #same substitution
#Fourth Normal Mode - same process as above
arg4 = sym.solve(KwM.subs(ω**2,ω_norm4)*A_matrix,(a1,a2,a3,a4)) #same method of generating arguments
mode_norm4 = A_matrix.subs(arg4).subs(a4,1) #same substitution
mode_norm1, mode_norm2, mode_norm3, mode_norm4 #Displaying ALL results of all normal modes
The above result contains the normal modes for our system
#Animating 4 different normal modes
#Defining plot
fig, ax = plt.subplots(figsize=(10,10))
#Defining arrays for initial x positions AND normal modes from part c
x_init = np.array([-9, -3, 3, 9])
norm1 = np.array([mode_norm1[0], mode_norm1[1], mode_norm1[2], mode_norm1[3]])
norm2 = np.array([mode_norm2[0], mode_norm2[1], mode_norm2[2], mode_norm2[3]])
norm3 = np.array([mode_norm3[0], mode_norm3[1], mode_norm3[2], mode_norm3[3]])
norm4 = np.array([mode_norm4[0], mode_norm4[1], mode_norm4[2], mode_norm4[3]])
#Converting sympy functions to numpy functions of normal frequencies from part b
ω1_func = lambdify((k, m), ω_norm1, 'math')
ω1 = ω1_func(1,1) #plugging in k = 1 and m = 1 to solve for ω1
ω2_func = lambdify((k, m), ω_norm2, 'math')
ω2 = ω2_func(1,1) #plugging in k = 1 and m = 1 to solve for ω2
ω3_func = lambdify((k, m), ω_norm3, 'math')
ω3 = ω3_func(1,1) #plugging in k = 1 and m = 1 to solve for ω3
ω4_func = lambdify((k, m), ω_norm4, 'math')
ω4 = ω4_func(1,1) #plugging in k = 1 and m = 1 to solve for ω4
#Plotting each mass
m1_norm1, = ax.plot(x_init[0],1,'s',color='red', markersize = '25') #Plot a single point representing m1 of normal mode 1
m2_norm1, = ax.plot(x_init[1],1,'s',color='red', markersize = '25') #Plot a single point representing m2 of normal mode 1
m3_norm1, = ax.plot(x_init[2],1,'s',color='red', markersize = '25') #Plot a single point representing m2 of normal mode 1
m4_norm1, = ax.plot(x_init[3],1,'s',color='red', markersize = '25') #Plot a single point representing m3 of normal mode 1
m1_norm2, = ax.plot(x_init[0],2,'s',color='blue', markersize = '25') #Plot a single point representing m1 of normal mode 2
m2_norm2, = ax.plot(x_init[1],2,'s',color='blue', markersize = '25') #Plot a single point representing m2 of normal mode 2
m3_norm2, = ax.plot(x_init[2],2,'s',color='blue', markersize = '25') #Plot a single point representing m2 of normal mode 2
m4_norm2, = ax.plot(x_init[3],2,'s',color='blue', markersize = '25') #Plot a single point representing m3 of normal mode 2
m1_norm3, = ax.plot(x_init[0],3,'s',color='green', markersize = '25') #Plot a single point representing m1 of normal mode 3
m2_norm3, = ax.plot(x_init[1],3,'s',color='green', markersize = '25') #Plot a single point representing m2 of normal mode 3
m3_norm3, = ax.plot(x_init[2],3,'s',color='green', markersize = '25') #Plot a single point representing m2 of normal mode 3
m4_norm3, = ax.plot(x_init[3],3,'s',color='green', markersize = '25') #Plot a single point representing m3 of normal mode 3
m1_norm4, = ax.plot(x_init[0],4,'s',color='purple', markersize = '25') #Plot a single point representing m1 of normal mode 4
m2_norm4, = ax.plot(x_init[1],4,'s',color='purple', markersize = '25') #Plot a single point representing m2 of normal mode 4
m3_norm4, = ax.plot(x_init[2],4,'s',color='purple', markersize = '25') #Plot a single point representing m2 of normal mode 4
m4_norm4, = ax.plot(x_init[3],4,'s',color='purple', markersize = '25') #Plot a single point representing m3 of normal mode 4
#Defining animation function for time t
def animate(t):
#Calculating positions
x_norm1 = x_init + norm1*(np.cos(ω1*t)+np.sin(ω1*t)) #positions for all masses of normal mode 1
x_norm2 = x_init + norm2*(np.cos(ω2*t)+np.sin(ω2*t)) #positions for all masses of normal mode 2
x_norm3 = x_init + norm3*(np.cos(ω3*t)+np.sin(ω3*t)) #positions for all masses of normal mode 3
x_norm4 = x_init + norm4*(np.cos(ω4*t)+np.sin(ω4*t)) #positions for all masses of normal mode 4
#Updating positions
#Normal Mode 1
m1_norm1.set_xdata(x_norm1[0])
m2_norm1.set_xdata(x_norm1[1])
m3_norm1.set_xdata(x_norm1[2])
m4_norm1.set_xdata(x_norm1[3])
#Normal Mode 2
m1_norm2.set_xdata(x_norm2[0])
m2_norm2.set_xdata(x_norm2[1])
m3_norm2.set_xdata(x_norm2[2])
m4_norm2.set_xdata(x_norm2[3])
#Normal Mode 3
m1_norm3.set_xdata(x_norm3[0])
m2_norm3.set_xdata(x_norm3[1])
m3_norm3.set_xdata(x_norm3[2])
m4_norm3.set_xdata(x_norm3[3])
#Normal Mode 4
m1_norm4.set_xdata(x_norm4[0])
m2_norm4.set_xdata(x_norm4[1])
m3_norm4.set_xdata(x_norm4[2])
m4_norm4.set_xdata(x_norm4[3])
#Returning the motion of all 16 masses
return m1_norm1, m2_norm1, m3_norm1, m4_norm1, m1_norm2, m2_norm2, m3_norm2, m4_norm2, m1_norm3, m2_norm3, m3_norm3, m4_norm3, m1_norm4, m2_norm4, m3_norm4, m4_norm4
#Creating list of time values in seconds
tlist = np.linspace(0,20,200)
#Running animation function
ani = FuncAnimation(fig, animate, tlist, interval=50, blit=True)
#Plotting and animating!
plt.title('Four Normal Oscillation Modes of Four Equal Masses Connected by Equal Springs')
plt.xlabel('X Position (m)')
plt.xlim(-12,12)
plt.ylim(0,5)
#Giving new tick marker values for y-axis - code modified from https://stackoverflow.com/questions/11244514/modify-tick-label-text
fig.canvas.draw()
labels = [item.get_text() for item in ax.get_yticklabels()]
labels[0] = ''
labels[1] = 'Normal Mode 1'
labels[2] = 'Normal Mode 2'
labels[3] = 'Normal Mode 3'
labels[4] = 'Normal Mode 4'
labels[5] = ''
ax.set_yticklabels(labels)
#Closing plot and converting HTML
plt.close()
HTML(ani.to_jshtml())